
#include "../../include/matrixFreePDE.h"

template <int dim, int degree>
void  MatrixFreePDE<dim,degree>::computeIntegral(double& integratedField, int index, std::vector<vectorType*> variableSet) {
	QGauss<dim>  quadrature_formula(degree+1);
	FE_Q<dim> FE (QGaussLobatto<1>(degree+1));
	FEValues<dim> fe_values (FE, quadrature_formula, update_values | update_JxW_values | update_quadrature_points);
	const unsigned int   n_q_points    = quadrature_formula.size();
	std::vector<double> cVal(n_q_points);

	// constraintsDirichletSet[index]->distribute(*variableSet[index]);
	// constraintsOtherSet[index]->distribute(*variableSet[index]);
	// variableSet[index]->update_ghost_values();

	typename DoFHandler<dim>::active_cell_iterator cell= this->dofHandlersSet[0]->begin_active(), endc = this->dofHandlersSet[0]->end();

	double value = 0.0;

	for (; cell!=endc; ++cell) {
		if (cell->is_locally_owned()){
			fe_values.reinit (cell);

			fe_values.get_function_values(*variableSet[index], cVal);

			for (unsigned int q=0; q<n_q_points; ++q){
				value+=(cVal[q])*fe_values.JxW(q);
			}
		}
	}

	value=Utilities::MPI::sum(value, MPI_COMM_WORLD);

	//   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0){
	//   std::cout<<"Integrated field: "<<value<<std::endl;
	//   }

	integratedField = value;
}

//-----------------------------------

//update RHS of each field
template <int dim, int degree>
void MatrixFreePDE<dim,degree>::computeIntegralMF(double& integratedField, int index, const std::vector<vectorType*> variableSet){
  //log time
  computing_timer.enter_section("matrixFreePDE: computeIntegralMF");

  integrated_var = 0.0;
  integral_index = index;

  //call to integrate and assemble
  matrixFreeObject.cell_loop (&MatrixFreePDE<dim,degree>::getIntegralMF, this, residualSet, variableSet);

  integratedField=Utilities::MPI::sum(integrated_var, MPI_COMM_WORLD);

  //end log
  computing_timer.exit_section("matrixFreePDE: computeIntegralMF");
}

template <int dim, int degree>
void MatrixFreePDE<dim,degree>::getIntegralMF(const MatrixFree<dim,double> &data,
                                        std::vector<vectorType*> &dst,
                                        const std::vector<vectorType*> &src,
                                        const std::pair<unsigned int,unsigned int> &cell_range) {

	dealii::FEEvaluation<dim,degree,degree+1,1,double> var(data,0);

    //loop over cells
    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell){
		var.reinit(cell);
		var.read_dof_values_plain(*src[0]);
		var.evaluate(true, false, false);

        unsigned int num_q_points = var.n_q_points;

		dealii::AlignedVector<dealii::VectorizedArray<double> > JxW(num_q_points);
		var.fill_JxW_values(JxW);

        //loop over quadrature points
        for (unsigned int q=0; q<num_q_points; ++q){


			dealii::VectorizedArray<double> val = var.get_value(q);
			assembler_lock.acquire ();
			for (unsigned i=0; i<val.n_array_elements;i++){
				integrated_var += val[i]*JxW[q][i];
			}
			assembler_lock.release ();

        }


    }
}

#include "../../include/matrixFreePDE_template_instantiations.h"
